home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MINV.FOR < prev    next >
Text File  |  1985-11-29  |  4KB  |  169 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MINV
  5. C
  6. C        PURPOSE
  7. C           INVERT A MATRIX
  8. C
  9. C        USAGE
  10. C           CALL MINV(A,N,D,L,M)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A - INPUT MATRIX, DESTROYED IN COMPUTATION AND REPLACED BY
  14. C               RESULTANT INVERSE.
  15. C           N - ORDER OF MATRIX A
  16. C           D - RESULTANT DETERMINANT
  17. C           L - WORK VECTOR OF LENGTH N
  18. C           M - WORK VECTOR OF LENGTH N
  19. C
  20. C        REMARKS
  21. C           MATRIX A MUST BE A GENERAL MATRIX
  22. C
  23. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  24. C           NONE
  25. C
  26. C        METHOD
  27. C           THE STANDARD GAUSS-JORDAN METHOD IS USED. THE DETERMINANT
  28. C           IS ALSO CALCULATED. A DETERMINANT OF ZERO INDICATES THAT
  29. C           THE MATRIX IS SINGULAR.
  30. C
  31. C     ..................................................................
  32. C
  33.       SUBROUTINE MINV(A,N,D,L,M)
  34.       DIMENSION A(1),L(1),M(1)
  35. C
  36. C        ...............................................................
  37. C
  38. C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
  39. C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
  40. C        STATEMENT WHICH FOLLOWS.
  41. C
  42. C     DOUBLE PRECISION A,D,BIGA,HOLD
  43. C
  44. C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
  45. C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
  46. C        ROUTINE.
  47. C
  48. C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
  49. C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  ABS IN STATEMENT
  50. C        10 MUST BE CHANGED TO DABS.
  51. C
  52. C        ...............................................................
  53. C
  54. C        SEARCH FOR LARGEST ELEMENT
  55. C
  56.       D=1.0
  57.       NK=-N
  58.       DO 80 K=1,N
  59.       NK=NK+N
  60.       L(K)=K
  61.       M(K)=K
  62.       KK=NK+K
  63.       BIGA=A(KK)
  64.       DO 20 J=K,N
  65.       IZ=N*(J-1)
  66.       DO 20 I=K,N
  67.       IJ=IZ+I
  68.    10 IF( ABS(BIGA)- ABS(A(IJ))) 15,20,20
  69.    15 BIGA=A(IJ)
  70.       L(K)=I
  71.       M(K)=J
  72.    20 CONTINUE
  73. C
  74. C        INTERCHANGE ROWS
  75. C
  76.       J=L(K)
  77.       IF(J-K) 35,35,25
  78.    25 KI=K-N
  79.       DO 30 I=1,N
  80.       KI=KI+N
  81.       HOLD=-A(KI)
  82.       JI=KI-K+J
  83.       A(KI)=A(JI)
  84.    30 A(JI) =HOLD
  85. C
  86. C        INTERCHANGE COLUMNS
  87. C
  88.    35 I=M(K)
  89.       IF(I-K) 45,45,38
  90.    38 JP=N*(I-1)
  91.       DO 40 J=1,N
  92.       JK=NK+J
  93.       JI=JP+J
  94.       HOLD=-A(JK)
  95.       A(JK)=A(JI)
  96.    40 A(JI) =HOLD
  97. C
  98. C        DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS
  99. C        CONTAINED IN BIGA)
  100. C
  101.    45 IF(BIGA) 48,46,48
  102.    46 D=0.0
  103.       RETURN
  104.    48 DO 55 I=1,N
  105.       IF(I-K) 50,55,50
  106.    50 IK=NK+I
  107.       A(IK)=A(IK)/(-BIGA)
  108.    55 CONTINUE
  109. C
  110. C        REDUCE MATRIX
  111. C
  112.       DO 65 I=1,N
  113.       IK=NK+I
  114.       HOLD=A(IK)
  115.       IJ=I-N
  116.       DO 65 J=1,N
  117.       IJ=IJ+N
  118.       IF(I-K) 60,65,60
  119.    60 IF(J-K) 62,65,62
  120.    62 KJ=IJ-I+K
  121.       A(IJ)=HOLD*A(KJ)+A(IJ)
  122.    65 CONTINUE
  123. C
  124. C        DIVIDE ROW BY PIVOT
  125. C
  126.       KJ=K-N
  127.       DO 75 J=1,N
  128.       KJ=KJ+N
  129.       IF(J-K) 70,75,70
  130.    70 A(KJ)=A(KJ)/BIGA
  131.    75 CONTINUE
  132. C
  133. C        PRODUCT OF PIVOTS
  134. C
  135.       D=D*BIGA
  136. C
  137. C        REPLACE PIVOT BY RECIPROCAL
  138. C
  139.       A(KK)=1.0/BIGA
  140.    80 CONTINUE
  141. C
  142. C        FINAL ROW AND COLUMN INTERCHANGE
  143. C
  144.       K=N
  145.   100 K=(K-1)
  146.       IF(K) 150,150,105
  147.   105 I=L(K)
  148.       IF(I-K) 120,120,108
  149.   108 JQ=N*(K-1)
  150.       JR=N*(I-1)
  151.       DO 110 J=1,N
  152.       JK=JQ+J
  153.       HOLD=A(JK)
  154.       JI=JR+J
  155.       A(JK)=-A(JI)
  156.   110 A(JI) =HOLD
  157.   120 J=M(K)
  158.       IF(J-K) 100,100,125
  159.   125 KI=K-N
  160.       DO 130 I=1,N
  161.       KI=KI+N
  162.       HOLD=A(KI)
  163.       JI=KI-K+J
  164.       A(KI)=-A(JI)
  165.   130 A(JI) =HOLD
  166.       GO TO 100
  167.   150 RETURN
  168.       END
  169.